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ABSTRACT 

Based on a basic principle of orbital 
resonance, we have identified a huge deficit of 
solar radiation induced by the combined 
amplitude and frequency modulation of the 
J Earth's obliquity as possibly the causal 
| mechanism for ice age glaciation. Including this 
j modulation effect on solar radiation, we have 
performed model simulations of climate change 
for the past 2 million years. Simulation results 
show that: (1) For the past 1 million years, 

: temperature fluctuation cycles were dominated 

^ by a 100-Kyr period due to amplitude-frequency 
resonance effect of the obliquity. (2) From 2 to 
1 million years ago, the amplitude-frequency 
1 interactions of the obliquity were so weak that 
they were not able to stimulate a resonance 
effect on solar radiation. (3) Amplitude and 
frequency modulation analysis on solar radiation 
provides a series of resonance in the incoming 
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solar radiation which may shift the glaciation 
cycles from 41-Kyr to 100-Kyr about 0.9 million 
years ago. These results are in good agreement 
with the marine and continental paleochmate 
records; Thus, the proposed climate response to 
the combined amplitude and frequency 
modulation of the Earth's obliquity may be the 
key to understanding the glaciation puzzles in 
paleoclimatology. 


L INTRODUCTION 

Incoming solar radiation, generally called 
"insolation,” is a key factor in climate change 
studies. Deviations of the incoming solar 
radiation depend upon astronomical parameters 
which characterize the orbit of the Earth around 
the Sun and the Earth's axis of rotation. An 
equation for the insolation deviations was 
formulated and solved by Milankovitch [1] and 
Vemekar [2] using numerical procedure. In the 
conventional insolation theory for climate 
change, the Milankovitch insolation equation is 
considered as an external forcing function in 
climate models [3]. The periodicities of the 
Milankovitch insolation cycles always depend 
upon three orbital parameters: (1) orbital 

eccentricity for the 100-Kyr cycle, (2) 
obliquity e for the 41-Kyr cycle, and (3) 
precession index e sin co for the 23-or 19-Kyr 
cycle (a> is the argument of perihelion). 

For several decades, geophysicists have 
been trying to develop physical mechanisms in 
which the 100-Kyr cycles in the ellipticity of the 
Earth's orbit around the Sun shift the pattern of 
insolation, triggering the growth and decay of 
great ice sheets in the high latitudes of the 
Northern Hemisphere. This effort has 
encountered great difficulty: The 100-Kyr 
orbital eccentricity cycle is too small in 
amplitude and too late in phase [3, 4] to produce 


the 100-Kyr glaciation cycle during the last 1 
million years. This is one of the most perlexing 
and enduring puzzles in paleoclimatology. 

The second puzzle is how to explain the 
ice sheet cycles in the Northern Hemisphere 
from 1 million to 2 million years ago when the 
global ice volume and deep sea temperature 
varied at an almost metronomic 41-Kyr period 
of the Earth's obliquity. Thirdly, we still do not 
understand why the 100-Kyr ice age cycle 
become dominant about 0.9 million years ago. 
Numerous hypotheses and models have been 
proposed to explain these glaciation puzzles. 
Proposed explanations for these puzzles have 
invoked: (1) stochastic resonance of 

eccentricity forcing [5], (2) internal oscillations 
of the climate system near the 100-Kyr period 
that can get phase-locked to external 
(eccentricity) forcing [6], (3) high nonlinear 
response of climate system to weak forcing by 
the orbital eccentricity [3, 7-13], (4) variations 
in the inclination of the Earth’s orbital motion 
[14-16] and (5) climate system with three steady 
states and a set of predefined rules for moving 
between them [17]. However, the physical 
mechanisms of the glaciation cycles are 
unknown. 

From dynamics point of view, Liu [1 8,19] 
has shown that the 100-Kyr periodicity of the 
ice ages is not due to eccentricity at all, but to 
frequency variations in the obliquity, which is 
the tilt angle of the Earth’s rotation axis from its 
orbital plane. Indeed, Rubincam [20,21] has 
obtained an analytical insolation equation which 
shows that neither the main pacemaker of the 
ice age e, nor Milankovitch precession index e 
sin a> appear as terms in the equation, but 
obliquity e does. Rubincam's analytical 
insolation equation is developed in terms of 
spherical harmonies. It is a convenient formula 
for computing insolation forcing when using 
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numerical climate models. In a series of papers 
[22-25], we have used numerical climate models 
to demonstrate that frequency modulation of the 
Earth's obliquity is responsible and accountable 
for major climate changes during the past 2 
million years. For a physical mechanism to be 
considered appropriate and convincing, we have 
computed variations in the solar energy flux at 
the top of the atmosphere and identified that the 
insolation flux deficit is the physics behind the 
ice age glaciation. 

We should emphasize that our interest in 
this study does not lie in the mathematical 
development of celestial mechanics ■ and 
sophisticated climate modeling. Rather, we are 
concerned with the necessary, fundamental 
concepts in wave physics to understand the 
climate cycles as seen in oxygen isotopic $ 8 0 
measurements from ocean sediment cores [26- 
29] and biogenic silica measurements from 
continental sediment cores in Lake Baikar [30]. 

Finally, we note that the obliquity 
modulation theory for climate change in this 
study is developed in order to relate the 
variations of the Earth's obliquity and its derived 
frequency to climate. However, changes in the 
solar constant due to natural evolution of the 
Sun as a star or related to sunspots will not be 
part of the theory. Different astronomical 
theories of paleoclimate, such as those of 
LeVerrier [31], Croll [32], Milankovitch [1], 
Vemekar [2] Berger and Loutre [33] and 
Rubincam [20,21], fall under such a definition. 


2. OBLIQUITY VARIATIONS 

Lunar and solar gravitational torques 
acting on the equatorial bulge of the Earth cause 
its spin axis to precess about the instantaneous 
orbital normal. The lunar and solar torques 


together produce a precession of the spin axis of 
the Earth at a rate of 50.38 arcsec yr l [34,35]. 
Once the present spin axis direction is known 
and orbital element histories are given, the 
obliquity history can be constructed by applying 
the methods, developed by Brouwer and van 
Woerkom [36]. In this way the variations in the 
Earth’s obliquity have been determined by [2]: 

e(t) = e* cos^ + K)t +64 + a] 

i-i 

-X c s N i cos(2[(s ( + K )t + 5 j +“]} 

- XX ' CijNiNi cosf( Sj + Sj + 2 K)t + 6 , + 6 , + 2 a] 

i-I j>1 

-XXG.N.N, c<J( Si - Sj )t + 8, - S] 

i-l j>l 

In equation (1), t is the time and all constants 
were given in [2]. The variations in the 
obliquity were computed from equation (1) foi 
the past 4000 Kyr with an interval of 1 Kyr b> 
Vemekar [2]. His results showed that the 
obliquity of the Earth varied in the past 4000 
Kyr only from 24.51° to 22.10°, with an average 
period of 41 Kyr. The obliquity period of 41 
Kyr has long been recognized as one for the 
Milankovitch glaciation cycles as seen in 
geological climate records. However, equation 
(1) represents a waveform with variable phase 
angles. Phase modulation may induce 
asymmetry of the obliquity cycles, permitting 
transient or time-dependent unstable climate 
features to be detected. This seems to be the 
pitfall in the Milankovitch theory for climate 
change by which we were led into interpretation 
difficulties of climate response to orbital forcing 
[3,37-40]. Probably the most important feature 
through which the orbital imprint may be 
unambiguously recognized in ancient geological 
records is the phase modulation of the obliquity 
component. Milankovitch theory is virtually 
valueless to detect this property in geological 
records. Actually, investigation of the phase 
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modulation effect of the obliquity on climate 
change requires complex representation of 
equation (1). 


<(>(t) = tan 1 


^AiSin^it + ^i) 



XA i co^y,t + ^) 
1*1 


(5) 


3. FREQUENCY MODULATION OF THE 
OBLIQUITY 

In order to compute the frequency 
modulation of the obliquity, the complex 
representation of e(t) can be expressed by a 
complex waveform 


Liu, [18] has suggested that climate 
response to obliquity forcing may be determined 
by the frequency variation of the obliquity. In 
this case the instantaneous frequency rai(t) of the 
phase-modulated obliquity is defined as the 
derivative of the phase: 


Ae(t) = e(t) - e" 

= 2X cos (yi t+ ^) (2) 

i-l 

= r(t)exp(V-18(t)) 


where A if t- x and Q are the amplitude, mean 
frequency and phase for each component, 
respectively. The initial value of s is 
23.230556° and t = 0 refers to AD 1950. All 

constants for i = 1,2,3* were given in 

[41,42]. The complex vector in equation (2) has 
a magnitude r(t) and angle 0(t) where r(t) is 
slowly varying and 0(t) varies rapidly. By 
setting the phase function 

0(t) = TD c t + $(t) (3) 

where w c is the mean obliquity frequency and' 
(j)(t) is the phase angle of the obliquity, then 


1 1 r Ae(x) 

Mt) = A e (t)®- = -L7r7 dT 


1 r» Ae(t - 1 ) 


dx 


, 7C 

= XA i sin(y i t + i; i ) 


(4) 


where Ag(t) is the convolution of Ae(t) and 
l/(7it). Therefore, 


where ra c = 31.079137 (arc sec yr 1 ) and d<t»/dt is 
the frequency deviation. Denoting 


Ata(t) = 


W) 

dt 


(7) 


we may define a phase factor p(t) as the ratio of 
the change rate of the phase function 0 of the 
obliquity to the mean frequency of the obliquity. 
Thus 


l d0 . Atu(t) 

^ (t) = ^7¥ =1+ Wc 


1 d 

. -i 

XAjSin^t + t;,) 

i-J 

tn c dr 

^ tan 

^cos^t^) 

_ i-l -L 


which regulates climate model response to the 
incoming solar radiation. 

The scientific basis of this study is the 
time series of the instantaneous frequency of the 
obliquity t3i(t) in equation (6). Variations of 
t 3 i(t) in time with their correspondence period 
variations are shown in Figure 1. Figure 1 
shows that the period of the obliquity varied 
from about 37 Kyr to 44 Kyr for the past 1 
million years. This lengthening/shortening 
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Figure 1. Frequency and period variations of the 
Earth’s obliquity. 



Figure 2. Spectral amplitude of the obliquity 
frequency variations. 


effect of the obliquity period is closely 
associated with major paleoclimate events [18]. 


indicates that 185-Kyr climate cycles could be 
recorded in geological data, as suggested by 
Beaufort [40]; (2) the deep split of the 100-Kyr 
peak into two peaks at about 120-Kyr and 88- 
Kyr is of particular interest because the irregular 
ice age cycles varied from about 80 Kyr to 120 
Kyr over the past 1 million years [19]; (3) the 
major peaks are seen at 100-Kyr and 185-Kyr, 
and the peaks at 41 -Kyr and 26- to 23-Kyr have 
relatively low intensities with minor 
significance; and (4) the peak at 45-Kyr was 
recorded by the core data of VI 9-28 and VI 9- 
29 in the Pacific Ocean. 

4. A TEST MODEL: FREQUENCY 

MODULATION EFFECT OF THE 
OBLIQUITY ON CLIMATE CHANGE 

The hypothesis of frequency forcing of 
the obliquity for climate change is that glacial 
and interglacial intensities are related to the 
frequency variation of the obliquity. Ice ages 
started after the instantaneous frequency tDj(t) 
reached a low value. The physical concept was 
described by Liu [19], but its validity needs to 
be tested by simple but realistic climate models. 
A class of zero-dimensional energy balance 
models (EBM) for climate change with surface 
temperature T as the main variable was 
developed as one of the fundamental climate 
models [5, 43-46]. For a simplified version of 
EBM, the equation of the model is 

c^ = Q (0[l-a(T)] + Em (9) 


Also, Liu [18] has reported a prominent 
100-Kyr peak in the spectrum of the xd\ (t) time 
series which may be related to the. anomalous 
100-Kyr peak in the spectrum. An 

improved spectral amplitude of TUj(t) is shown in 
Figure 2. The physical interpretation of this 
spectrum may be made as follow^: (1) it 


where c is the thermal capacity of the Earth, 
Q(t) is the incoming solar radiation, a(T) is the 
globally averaged albedo, and E(T) is the 
outgoing surface radiation. The formula for the 
albedo is [5]. 

a(T) = d 1 -d 2 tanh[d J (T-T.)] (10) 
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where di, d 2 and d 3 are constants. E (T) can be 
parameterized as a linear function of T through 
two empirical parameters A and B such that 

E(T) = ~(A +■ BT) (11) 

Sojar intensity Q(t) reaching the Earth may be 
modulated by the frequency variability of the 
obliquity [18,19], which is a function of the rate 
of change of the phase function 0 as defined by 
equat'pn (8). Thus 

Q(t) = n(t)Q. (12) 

where Q 0 is one quarter of the solar radiation 
intensity. The parameter fj(t) in equation (12) 
allows us to introduce explicit variations in the 
solar input. Accordingly, the equation of the 
EBM becomes 


dT 1 d 


i«»1 


C dt w c dt 

tan 

Z A i cos (yi t+ 0 

N i-1 ' 



• {l - d, + d 2 tanh[d 3 (T - T.)]| - A - BT 


The constants in equations (10) and (11) for 
albedo and the outgoing surface radiation are 
chosen to yield roughly the limits of the 
temperature range exhibited by the glacial- 
interglacial oscillations of the Pleistocene and, 
therefore, they could be interpreted as 
representing the glacial and interglacial states of 
the Earth's climate system. Using other 
parameterizations for the albedo and outgoing 
surface radiation [43, 44, 47], the model is able 
to produce slightly different results for 
comparison. 

The estimates of the geophysical 
constants in equation(13) were given in [5] and 
all numerical values of orbital constants were 
given in [41, 42] for i — 1,2,3, We have 
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computed temperature T as a function of time 
from equation (13). With the above choices of 
the orbital constants and geophysical constants, 
we have obtained that the range of variations in 
T is about 11° C in paleoclimatic variations 
assuming a linear relationship between $ 8 0 and 
the surface temperature [48,49]. However, the 
8°C temperature change commonly inferred 
from core isotopic records [50-52] is relatively 
smaller than the recent estimates of the 15° C 
temperature change from central Greenland 
[53]. 

The time-dependent solution of equation 
(13) are derived from the lengthening/shortening 
of the obliquity cyclic period. In order to 
display the relationship between the temperature 
T and fi! 8 0 variations, we have determined the 
time-dependent changes in the amplitude of the 
100-Kyr frequency component in both the T(t) 
and the S? 8 0 time series. The method of 
complex demodulation [54] allows us to make 
such a determination. Complex demodulation is 
carried out by first calculating the product of T (l) 
exp (i 2 k ft) where i is the square root of -1 and 
f is the 100-Kyr frequency at which we wish to 
determine the time-dependent changes in 
amplitude. Complex demodulation of the 
simulated temperature T(t) and the complex 
demodulation of the $ s 0 times series [55] are 
shown in Figure 3, indicating the general 
similarity between the results of simulaions and 
the observation data. Therefore the test model 
indicates that equation (13) can be used to 
generate a 100-Kyr modulating cycle in 
insolation. Now, our critical problem is to 
explain, physically and mathematically, the 
process through which equation (13) will 
modulate the insolation as a causal mechanism 
for climate change. 

Physically, obliquity variation equations 
(1) and (2) represent an oscillaton signal whose 
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Figure 3. The solid curve denotes complex 
demodulation of the simulated temperature 
fluctuations caused by frequency variations in the 
obliquity. The dashed curve denotes complex 
demodulation of the 5 I8 0 time series from [29]. 


amplitude is slowly varying with respect to 
frequency. Therefore application of obliquity 
variation as a component of insolation forcing 
for climate change should involve both the 
amplitude and frequency modulation, and their 
possible interactions or resonances. Frequency 
modulation is fully described by equations (4), 
(6) and (8). Since these frequency variation 
equations can be used to determine the 
lengthening/shortening of the cyclic period of 
the obliquity which are related to the prolonged 
glaciation and rapid melting of ice age ice sheets 
in the northern hemisphere, it seems justifiable 
to assume that climate model response to the 
incoming solar radiation is modulated by 
frequency variation of the obliquity. 

Another question that might be asked is: 
Why should frequency modulation be 
considered as a forcing in climate models? We 
consider frequency modulation of the oblquity 
as a component of insolation forcing for climate 
change for the following reasons: (1) Its 
resonances with the insolation could induce 


pulses in the incoming solar radiation; climate 
response to these insolation pulses would cause 
major climate changes. (2) The 
lengthening/shortening of the cyclic periods of 
the obliquity, determined by equation (8), could 
be observed. (3) It can be recognized by 
performing modulation analysis on the complex 
demodulation of the geological data. (4) Finally, 
the frequency-insolation resonance effect seems 
to disappear before the mid-Pleistocejie at about 
900 kyr B.P. Our theoretical justification of the 
assumptions in the model is gratifying because it 
brings these important characteristics to light. 

The motivation for equation (8) is the 
investigation of climate change via insolation. 
Ideally, equation (8) should provide (1) 
theoretical insights into the insolation and (2) a 
convenient formula for computing the insolation 
pulsation in the incoming solar radiation when 
using numerical climate models. The point of 
model simulation exercises is that equation (8) 
does seems to be suggestive of insolation 
modulation. Indeed, the test model indicates 
that equation (8) can be used to generate a 100- 
kyr modulating cycle in insolation. Therefore 
the simulation test of equation (8) as presented 
here is an example of the first motivation 
mentioned above. However, whether equation 
(8) fulfills the second motivation or not is yet to 
be established and justified. 

The 100-kyr result of the model 
simulation test of equation (8) does not mean 
that it is a separate forcing function for climate 
change. The purpose of equation (8) is to 
provide a mathematical expression of physical 
reasons for modification of the insolation as a 
forcing function. It should be noted that there is 
a deep-seated misconception about the 
Milankovitch insolation forcing for climate 
change. According to phase theory, application 
of the Milankovitch insolation equation as a 
forcing function for climate change is 




22 


Han-Shou Liu 


incomplete and inappropriate. Below it will be 
shown that frequency variation of the obliquity 
can induce insolation pulses and the pulse 
modulation of the insolation may provide an 
appropriate forcing mechanism to which the 
climate system would respond. 

5. INSOLATION EQUATIONS 

Milankovitch insolation function Q(t) is 
expressed in terms of orbital parameter 
variations. 

Q(t) = AR s Ae(t) + mA(esinco) (14) 

in which AR S is the difference of the insolation 
received at the top of the atmosphere over the 
caloric summer for 1° change in obliquity e(t), 
e sinco is the precession variation [2], 


where T y is the duration of the tropical years, S 
is the solar constant, and is the latitude. 
Variations in orbital parameters (e, e, esinco) are 
governed by 


Kyr, 41-Kyr and 23- or 19-Kyr, respectively. 
The peak-to-peak oscillation of the obliquity is 
about 2° [56], 

New mathematical analysis to improve 
the accuracy of the constant in equations (16), 
(17) and (18) have been developed in [2, 41, 57- 
62]. ON the basis of these unproved orbital 
constants, the time series from equation (14) for 
<jf=65°N is shown in Figure 4 [22]. Its wavelet 
spectrum shows that the Milankovitch- type 
insolation forcing provides an almost negligible 
contribution from eccentricity forcing for the 
100-kyr ice age cycles. A linear version of the 
Milankovitch insolation theory is inadequate to 
explain the 100-kyr glaciation cycle [3]. 

We propose to investigate how the 
climate system would respond to the amplitude- 
frequency resonance effect of the orbital 
parameters. In general, consideration of 
resonance in physics has proved to be enorously 
fruitful [63-69]. We argue here that it may 
provoke fundamental questions with regard to 
orbital forcing, insolation, and climate. In 
coming to grips with the physical meaning of the 
Milankovitch orbital forcing function, equation 
(14), as it relates to the nonlinear climate 


e = e. + 2^1 C0S (M + ^<) 

t-i 

(16) 

Mt) = Z A i co ^Yit + ^i) 

(17) 


i*i 


A(esina)) = ^Pj sir^ait + p,)- 0.0167 sin(l02') (18) 

It is noted that equation (17) is the derivation 
from a constant of the obliquity e° = 23.320556° 
and equation (18) is the derivation from the 
present-day value of esinco. It is well known 
that orbital eccentricity e, obliquity change As(t) 
, and precession index A(esinoo) oscillate at 100- 



Figure 4. Milankovitch insolation at 4> = 65°N 
and its wavelet spectrum. 
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problems, it is important to try to understand the 
concept of the instantaneous frequency of the 
insolation variations. 

The instantaneous frequency hypothesis is 
that glacial and deglacial intensities are related 
to the frequency variations of the insolation. Ice 
ages started after the instantaneous frequency of 
the insolation reached minimum values, and 
deglaciations occurred at times when the 
instantaneous frequency of the insolation 
reached maximum values. 

The time series from equation (14) is 
nonstationary and contains possible amplitude- 
frequency resonances in variations of the 
eccentricity, obliquity, and precession and their 
possible coupling interactions. Computing the 
instantaneous frequency of the insolation 
involves complex representation of the 
insolation function. The imaginary part of the 
function is actually the Hibert transform of the 
insolation Q(t) in equation (14). The Hibert 
transformation of insolation may provide a 
physical mechanism for climate change studies. 
However, we cannot make the Hibert transform 
for the Milankovitch insolation time series 
because the amplitude of the signal from 
equation (14) is not band-limited. It is possible 
that we can perform Hibert transform for the 
Rubincam insolation because the amplitude of 
the Rubincam insolation signal is in general 
band-limited. 

To a high degree of approximation, the 
insolation F s at any point on the Earth is given 
by [20] 

F s = F s -Hcos Vs [f) (19) 

where F s ° =1368 Wm' 2 is the average amount of 
sunlight impinging on the Earth, ro and r s are the 
average and instantaneous distance from the 


center of the Earth to the center of the Sun 
respectively, % is the solar zenith angel and 
H-l when cos W s >0 and H=0 when cos *Fs<0. 
When cos y/ s <0 the Sun is below the horizon 
and the flux is zero. 

Rubincam [20] has developed equation 

(19) in terms of ((p,X,rj,a,e,e,Q,to,M) where <p is 
latitude, X longitude, a the semimajor axis of the 
orbit, e the obliquity, Q the nodal position, co the 
argument of perihelion and M the mean 
anomaly, all these orbital parameters are 
referred to an Earth-fixed frame. The hour angle 
tj , measured from a fixed point, describes the 
Earth's rotation. 

In expressing insolation in terms of the 
Earth's orbital parameters, equations (19) 
becomes » 

p-0q-<° l M-modd 

•[(/-2/?)ra + (/ -2p + q}M+n{Q.-r[-'k^ 

which is the fundamental equation of insolation. 
In equation (20), the P[ m (sin <p) are the 
associated Legendre polynomials, F\ mp (s) are 
the inclination functions commonly used in 
celestial mechanics, Wi. 2 P , q (e) are eccentricity 
functions of the insolation and d\ are the 
coefficients of the insolation. 

Averaging over the diurnal and annual 
cycles, leaves only the terms with m = 0 and 
l-2p+q = 0. Furthermore, W n> .„ = 0 for all n > 
0 , so that only W 0t o = 1 + (l/2)e 2 +(3/8)e 4 
is non-zero. Finally, we note that di vanish for 
all odd / and decrease rapidly for all even /. 
Therefore, the truncation of l > 2 for equation 

(20) is a reasonable approximation for our 
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purpose. As a result, only the terms with / = 0, 
p = 0 and / = 2, p = 1, remain, and the zonal 
insolation to zero order in e is given by 


x {i + ^ A(sin<p) 

jsin 3 s(0-|| 


(21) 


where P 2 (situp) is the ordinary Legendre 
polynomial of degree 2, and terms in e 2 and e 4 
are dropped due to their small values compared 
to 1 . It should be noted that there is no e sin co 
terms in equation (21). Such terms depend on 
season, and the annual average is zero. But 
obliquity s does appear. By writing e = e(t) = e° 
+ As(t), dependence of the insolation on the 
obliquity Ae(t) becomes obvious. Combination 
of amplitude and frequency modulation of 
insolation may induce pulses in insolation. 


6. PULSATION OF THE INSOLATION 

Variations of the Rubincam insolation are 
computed from equations (17) and (21) for the 
past million years with an interval of 100 years. 
The results are shown in Figure 5. The real part 
of the complex vector of the insolation deviation 
A F/ (t) = F s *(t) - 222.125 from equation (21) 
can be expressed by 

AF* (t) = r(t)cos[e(t)] = r(t)cos[iD c (t) + £(t)] (22) 


where r(t) is the amplitude, 0(t) is the argument, 
w c is the mean frequency of the insolation 
variations, and q(t) is the phase angle. By 
convocation theory [70, p. 16] the imaginary 
part of insolation is defined by the Hibert 
transform 





Figure 5. Variations in Rubincam’s insolation. 



Urn* [in Kyr) 



Figure 6. Frequency variations in Rubincam’s 
insolation. 


Therefore the time derivative of the phase angle 
g is identical to the frequency variation of the 
insolation 


dg(t) 

dt 


d } 

—i arc tan 

'(AF*(t))| 

dt 1 

AF§ (t) j 


(24) 


The results of the phase modulation dg(t) / dt of 
Rubincam's insolation and their spectrum are 
shown in Figure 6. Figure 6 shows that 
frequency of the insolation varies with 100-kyr 
and 200-kyr timescales. 

Figure 6 may provide a physical 


t-x 
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mechanism of the obliquity to cross the 
threshold of insolation. The threshold model 
presented here is actually derived from a signal- 
to-noise model. A threshold system 
[17,22,69,71] for insolation is shown in Figure 
7. The insolation resonance consists of a 
threshold (shown by two horizontal dashed 
lines) and a subthreshold signal of the insolation 
at 65°N (thin curve with its frequency 
modulation signal (thick curve). The addition 
of the frequency to the insolation curve (thin) 
enables Rubincam's insolation curve to cross 
the threshold. Each time the insolation intensity 
plus its corresponding frequency (normalized) 
increase across the threshold, a pulse (positive 
or negative) can be written to the time series. 
The quantized output is a pulse modulation 
train, which _ reproduces more closely the 
characteristics of the analog input signal (thick), 
indicating the insolation resonance of orbital 
forcing as a paradigm for mode (phase or 
period) locking [72] in the climate system. 


To be more specific, the timing of 
insolation pulses essentially coincide^ with the 
occurrence of insolation resonances. These 
resonances with various magnitudes occurred at 
(-960, -910, -870, -800, -700, -630, -540, -500, 
-420, -330, : 220, -130, -10) kyr and (-930, 
-820, -760, -650, -600, -470, -390, -270, -187, 
-70) kyr for rapid ice sheet melting and deep 
glaciation, respectively (Figure 7b). 

The bipolar pulse modulation train c(t) in 
Figure 7b can be expressed by [70]. 

c(t) = E a , p ( t - t i)-Za i P(t-t j ) (25) 

i j 

where a-, and aj are the magnitude of each 
unequally spaced pulse, and P (t — t,) and P (t - 
tj) are the basic pulse shape in which t* and tj are 
the time indices of insolation resonance. 

The modulating function of the insolation 



Figure 7. Obliquity variations cause resonances and pulses in insolation: (a) 
Amplitude-frequency resonances in insolation, (b) Bipolar pulse modulation train of 
the insolation, (a* and aj are the magnitude of each pulse) 
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signal F s * ft) with a pulse modulation train eft) 
is Ff(t) eft). Therefore, from equations (17), 
(21), and (25) the structure of the pulse- 
modulated insolation can be expressed 

p[AF* (t)] = AF* (t) 1 + Za|P(t - t,) -2>jP(t - tj) 

(26) 

In trying to understand the behavior of 
insolation F/ ft) as a forcing function for 
climate change the analytical formula obtained 
from orbital consideration can be applied using 
pulse modulation analysis. This is what has 
been shown in equation (26). The time series 
of equation (26) is nonstationary and contains a 
series of amplitude-frequency resonances due to 
the Earth’s changing obliquity. 

The interpretation of the nonlinear 
forcing function (equation 26) is that the pulse 
modulation is actually the combined effect of 
physical processes that is absent from the 
climate models: (1) negative insolation pulses 
that induce ice ages when insolation is reduced 
during downward cycles, at times when the 
frequencies are lower, and (2) positive 
insolation pulses that induce rapid ice sheet 
melting when insolation is increased during 
upward cycles, at times when frequencies are 
higher. 

7. CLIMATE MODEL SIMULATIONS AND 
RESULTS 

(a) Energy Balance Models 

We have selected energy balance models 
to study the Baikal cycles and tf 7 8 0 variations 
because they are simple but realistic climate 
models. A class of zero-dimensional energy 
balance models (EBM) for climate change with 
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surface temperature T as the main variable was 
developed by Budyko [43], Sellers [44], 
Nicolis [45], North et al [46], and Matteucci 
[5] as one of the fundamental climate models. 
For a simplified version of EBM the equation of 
the model is 

c~ = p[aF* (t)|l - a(T)] + E(T) (27) 

where c is the thermal capacity of Earth, 
P[AF,*(t)] is the pulse-modulated insolation, 
a(T) is the globally averaged albedo, and Eft) is 
the outgoing surface radiation as defined by 
equations (10) and (1 1). 

Thus, from equations (26), (27), (10), and (11) 
we obtain 

= + - 0 _ ^ a * p ( t - *j) (28) 

- {l - d , + d 2 tanh[d 3 (T - T.)]} - A - BT 

The solution of equation (28) is a time 
series of T(t) as shown in Figure 8, 
superimposed on the time series of the biogenic 



Time (Kyr) 

Figure 8. Calculated temperature variations 
T(t) compared to the marine 5 l8 0 fluctuations 
and the continental biogenic silica cycles in 
Lake Baihal. 
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silica changes (Baikal Drilling Project (BDP)- 
96-2)) and tf 8 0 variations. The T(t) curve 
indeed shows large variations of the same order 
of magnitude as seen in paleociimate records 
[29, 30, 40]. To improve the accuracy of T(t) in 
amplitude and phase, the global model should 
be extended to one-dimensional formulation. 

For the 800 kyr time interval of overlap 
between the model results and proxy data, we 
have calculated their linear correlation and 
coefficients r. The Pearson coefficient r for the 
T(t) - tf 8 0 and T(t) - BDP correlation is about 
0.85 and 0.87, respectively. Thes.e are 
significant enough to demonstrate their linear 
correlation. However, a simple statistical 
correlation with the marine and terrestrial 
paleociimate records, although intriguing, is not 
a sufficient reason by itself to accept the 
hypothesis. To obtain a different view of the 
periodicity evolution, a wavelet spectrum 
analysis of the signals is needed [22]. We have 
performed wavelet spectral analysis of the time 
series of T(t) and BDP and found that the 
characteristic spectral features in the T(t) and 
BDP spectra are in good agreement with that of 
the £> 8 0 spectrum obtained by Bolton et al 
[73], Lau and Weng [72], and Liu and Chao 
[22]. It is noted that coherency analysis of T(t), 
tf 8 0, and BDP time series can also provide 
coherency spectra to understand their 
correlation in phase. 

(b) Ice sheet Models 

We have introduced the orbital resonance 
effect on insolation pulsation to examine the ice 
sheet climate models. The ice sheet climate 
models [19, 74-84] showed that the solar 
radiation variations due to orbital forcing are 
large enough to have produced ice ages. 
However, to make a small ice sheet grow to a 
critical large size or to make a critical large ice 


sheet decay rapidly, it is necessary that high 
rates of accumulation or ablation exist at least 
in the Northern Hemisphere during the 
glaciation and deglaciation phases. Therefore 
two predictions of ice sheet climate models 
should be tested against independent data: (I) 
higher rate of accumulation in the glaciation 
phase and (2) higher or rapid rate of ablation in 
the deglaciation phase. 

Three parameters ( a , b, and hj in 
Weertman's ice sheet model control the growth 
and decay of ice sheets: a is the accumulation 
rate, b is the ablation rate, and h s is the 
elevation of the snow line. The half width L of 
an ice sheet can be formulated In terms of a,b 
and h s , Nye [74, 75] and Weertman [76-78] 
showed that moderate changes in a, b, and h s 
can change L by large amounts and can even 
make it impossible for such an ice sheet to 
exist. They have estimated that ~ 100 years of 
greater than normal accumulation or ablation 
caused by weather fluctuations would be long 
enough to induce a small ice sheet to start 
growing to a critical size or to melt a large ice 
sheet rapidly. 

Under the low-low resonance conditions 
the negative insolation pulses last ~ 2 kyr 
longer than that in the normal glaciation cycles. 
Thus a small ice sheet would then expand into a 
critical large ice sheet as a result of increases in 
the accumulation rate for ~ 2 kyr. Rapid 
warming under high-high resonance conditions 
would increase the rate of change in T(t) and 
hence increase the rate of ablation [84]. The 
maximum insolation during rapid deglaciation 
interval shifted ~ 2 kyr in advance. Therefore 
higher ablation rates would set in and persist for 
~ 2 kyr longer than that in the normal 
deglaciation cycles. This seems to be the 
triggering mechanism for rapid melting of ice 
age ice sheets. 
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The center of the ice sheet on a continent 
in the Northern Hemisphere would shift 
southward as the ice sheet grows. The half 
width L of the ice sheet varies with time 
between the limits of about 0 < L < 1 000 km. 
On the northern edge of the ice sheet the snow 
line elevation is h( X -.jj(t), which varies in 
accordance with the variations in the ice sheet 
size. For reasonable variations in the values of 
a and b, Weertman’s model showed that h( X ^ 
i)(t) would fluctuate between the values of 0 
and 0.6 km if variations in the half width L are 
limited to the range 0-1 200km. Therefore the 
value of h( X ~.L)(t) would be zero at the 
occurrence of the positive insolation pulsation 
and 0.6 km at the occurrence, of the negative 
insolation pulsation. It is assumed in this simple 
way that the timing of the occurrence of the 
maximum and minimum values of the snow line 
elevation in the ice sheet is paced by the orbital 
resonance indices. Consequently, the elevation 
of snow line can be normalized and weighted 
by equation (26) [19]. 


The rate of change of the half width of ice sheet 
dL(t)/dt is given by [78]. 


dL(t) 2a 


dt 


1 + - 


3P 

4sh 


1/2 


(L(t) - ' 


1 + 


4s 

(*— l)( 0 — 2sL(t) 


(29) 


-1 


}T 2 (t) 


in which s is the snow line slope of the ice sheet 
and P - (4x)/(3pg), where t is the basal shear 
stress, p is the density of ice, and g is the 
gravitational acceleration. Frequency 

modulated insolation (equation (26)) will 
undoubtedily cause changes in a, b, h^ L )(t) in 
equation (29), as discussed above. To include 
the nonlinear processes in the ice sheet climate 
model induced by the resonances in insolation, 
we have modified and weighted the ice sheet 


parameters (a, b, h^L) (t) by equation (26), 
using the computation procedures developed by 
Liu [19]. Computer simulations were made 
incrementally by calculating AL(t)/At, where At 
is taken as 1 kyr. The values of P[AF s *(t)] 
were chosen for <j> = 65 °N. 

The results of computing equation (29) 
are shown in Figure 9. The time series of the 
half width of the ice sheet L(t) from equation 
(29) seems to mimic the time series 5 l8 0. 
Therefore introducing the frequency modulated 
insolation forcing function into the energy 
balance climate model and the ice sheet climate 
model yields reasonable matches to the 
observed climate records over the last 800 kyr 
(Figure 9). 

We have applied the simple but realistic 
energy balance and ice sheet climate models to 
demonstrate the amplitude-frequency resonance 
effect of the insolation on climate change at 
high latitudes. However, these climate models 
were derived decades ago and the climate 



Figure 9. Calculated ice sheet variations L(t) 
compared to the marine 6 18 0 fluctuations from 
Ocean Drilling Project Site 677. 
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changes more or less than the global average, 
depending mostly on latitude. While changes of 
~10°C would have occurred at high latitudes, 
changes of only -2°C would have occurred near 
the tropical latitudes. In order to uncover the 
detailed latitudinal signiture of the proposed 
physical mechanism it is necessary to apply the 
Ml general circulation models of climate. 


8. PHYSICAL MECHANISM OF 
OBLIQUITY FORCING 

Although the results from model 
calculations show that equation (26) may 
provide a key element for the interpretation of 
the terrestrial and marine climate change, the 
physical mechanism by which the climate 
system is forced is still not clear. If the climate 
system instead reacts to the regular insolation 
forcing and responds primarily to the obliquity's 
phase variations, then a physical mechanism 
should be recognizable. For a physical 
mechanism to be considered appropriate and 
convincing, we will need to compute the 
obliquity-forced insolation variations in terms of 
energy flux at the top of the atmosphere. 

When insolation is reduced during 
downward cycles of the obliquity, the total 
deficit of insolation on the top of the 
atmosphere in the Northern Hemisphere during 
a major (longer) glaciation period can be 
expressed by 

AE = I* AF* (t)dt = k(AE). (30) 

where tj and t 2 are the beginning and ending 
time of the longer glacial interval respectively. 
AF s *(t) = F/(t) - 222.125, A: is a numerical 
factor, and (AE)o is the average value of the 
total insolation deficit for a regular glacial 


interval. Computations of equation (30) for the 
past 10 major (longer) glacial time intervals 
yield 1.07 < k < 1.11. Therefore, under 
resonance conditions the total insolation deficit 
on the top of the atmosphere during the longer 
glacial intervals is -10% more than the average 
value (AE)o- This magnitude of insolation 
deficit seems sufficient to drive Earth into a 
deep ice age. Similarily, it can be shown that 
under resonance conditions the total insolation 
surplus on the top of the atmosphere in the 
Northern Hemisphere during the interglacial 
intervals is -10% more than that of the regular 
deglacial intervals. The results of the insolation 
intergral over time (equation (30)) are 
illustrated as shown in Figure 10: under the 
low-low resonance conditions, the negative 
insolation pulse (area bfgbf)) provides a 
physical mechanism of insolation forcing for 
prolonged glaciation, and under the high-high 
resonance conditions, the positive insolation 
pulse (area dfed’d) provides a physical 
mechanism of the insolation forcing for rapid 
deglaciation. Figure 10 shows that both 
negative and positive insolation pulses last -2 
kyr. It should be noted that there are two types 
of pulsation: intensity pulsation and duration 
pulsation. The negative and positive pulses in 
Figure 10 are defined as duration pulsation, not 
intensity pulsation. It is the duration pulsation 
of insolation that drives the climate to oscillate 
with 100-kyr cycle. 

Whatever the interpretation, the 
conclusion is that its realism can be measured 
by calculating the integral of equation (30) over 
time for having only an obliquity signal. This 
kind of insolation calculation can be used in the 
EBM to scale the resonance effect in insolation 
on climate change. 

Indeed, ignoring the pulsed source of 
insolation, F*(t) forces the model to give 0.4°C 
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Figure 10. Physical mechanism of insolation 
forcing for climate change for (a) low-low 
resonance conditions and (b) high-high 
resonance conditions. 

for the temperature range, which is in 
agreement with the results obtained by 
Matteucci [5]. However, modified by 
insolation pulses due to insolation resonances, 
equation (26) forces the model to give ~10°C 
for the temperature range. This shows that the 
magnitude and the formulation of the phase 
modulation influence of Rubincam's insolation 
for climate change are too big and important to 
be dismissed, even at a preliminary stage. 

As a forcing function, equation (26) 


shows that the climate system should respond 
not only to variations in the time-dependent 
insolation but also to a combined effect of 
changes in insolation and frequency of its quasi- 
sinusoidal variation. This is quite different from 
arguing that the system is highly nonlinear and 
threshold driven. From climatologist's point of 
view, interpretation of the forcing function is 
that the pulse modulation introduced here is 
actually a physical property of orbital dynamics 
that is absent from the climate models. 


9. WAVELET SPECTRUM OF THE 
MODELED TEMPERATURE 
VARIATIONS 

For the 800 kyr time interval of overlap 
between the model results and proxy data, we 
have calculated the Pearson's linear correlation 
coefficient r between the signals. The value of 
the Pearson's r is about 0.85, which is 
significant to support our conclusions. 
However, a simple statistical correlation with 
the paleoclimate record is not a sufficient 
reason by itself. For a new view of the 
periodicity _ evolution, a wavelet spectrum 
analysis of the signals is needed. [22]. 

Nonstational cycles abound in the time 
series T(t) from equation (28). The time- 
frequency spectral behavior of this model 
output data is of theoretical and practical 
interest. The most important feature of the T(t) 
time series is that its 100-kyr cycle is not stable 
through time. 

Advanced spectral tools such as 
Blackmann-Tukey, maximum entropy, and 
highly efficient Thomson technique have been 
applied to investigate the paleoclimatic 
variability using a set of four deep core data 
[85]. The dominant feature is the presence of a 
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period centered at 117 kyr. One of the best 
paleoclimate records is the $ 8 0 curve, obtained 
by Shackleton et al [29]. Recently, Bolton et 
al [73] and Lau and Weng [72] have performed 
wavelet spectral analysis on the S 8 0 time 
series and found that the 100-kyr cycle in the 
wavelet spectrum of the tf 8 0 data is time 
dependent. In this regard, we have performed 
wavelet spectral analysis on the T(t) time series 
to demonstrate that this model for Pleistocene 
climate provides a good fit to the time evolution 
characteristics in the wavelet spectrum of the 
observed tf 8 0 data. 

The wavelet transform of the time series 
T(t) with respect to the analyzing wavelet 
g(a,b)(0 can be defined as a convolution integral 
[ 86 ] 

Wg„. b) (t) = a- 1 ' 2 j>(t)g{^)dt (31) 

where a denotes the scale (dilation), b denotes 
the position (translation) and g° [(t - b)/a] is the 
complex conjugate of g (a ,b)(0 defined on the 
open real (a, b) half plane. The wavelet 
spectrum is thus displayed in the time- 
frequency domain, or the a-b space with the 
horizontal time axis b and the vertical frequency 
axis a. For our application we select the Morlet 
wavelet [87] which is a normalized, Gaussian- 
enveloped complex sinusoid with zero mean. 
We choose to examine the real part of the 
wavelet transform. It gives the amplitude 
undulation with the appropriate polarity and 
phase with respect to time owing to the 
symmetric nature of the real part of the Morlet 
wavelet as the kernel in integral equation (31). 
In contrast, the imaginary part, being 
antisymmetric in the kernel, gives the amplitude 
undulation as well, but imparts a 90° phase shift 
in time. In many application cases the modulus 
(combining real and imaginary parts) is 


preferred, but then the polarity phase 
information, which is important in our study, 
becomes absent. We use color contours such 
that amplitude peaks and troughs in horizontal 
successions in high contrast colors indicate the 
presence of strong oscillations ip the data, 
relative to the weaker and less significant 
background of low color contrast [22]. 

A wavelet spectrum of the time series 
T(t) thus developed is shown in Figure 1 1. One 
of the most important spectral features in Figure 
11 is that the 100-kyr cycle is not stable 
through time. After -650 kyr a near 100-kyr 
cycle dominates the spectrum, reflecting 
stronger pulses in the modulated radiation 
deviation. Prior to this time, in the region of 
time-frequency space near -850 kyr to -650 kyr, 
a near 80-kyr cycle emerges, in addition to a 
weakened 120-kyr cycle. Obviously, this is 
attributable to the strength and timing of the 
amplitude-frequency resonances of the 
obliquity. In this way, pulse modulation of the 
incoming solar radiation due to the amplitude- 
frequency resonance of the obliquity may cause 
the 100-kyr climate cycle. 

Indeed, the model results show a 
transition from dominant 41 -kyr to dominant 
100-kyr period about 900 kyr ago, as seen in 
the paleoclimatic record. Therefore the wavelet 
spectrum of the model’s results (Figure 11) 
explains clearly the periodicity evolution of the 
paleoclimatic fluctuations over the past 2 
million years. A simple Fourier analysis would 
suffice to prove that the paleoclimatic record is 
nonstationaiy over the whole Pleistocene. 
However, it can not demonstrate its periodicity 
evolution. 

The time-frequency spectral behavior of 
the model output as shown in Figure 1 1 is in 
good agreement with that of the tf 8 0 record 
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Figure 11. (a) Temperature variations and bipolar insolation pulses from the model 
and (b) wavelet spectrum of the temperature variations for the past 2 million years. 


[72,73]. We note that since the publication of 
the paper by Hays et al [88], numerous 
modeling efforts have attempted to explain the 
relation between Milankovitch insolation 
forcing and climate change. Most of these 
modeling studies have focused on the origin of 
the 100-kyr cycle. What these models do 
confirm is that climate response to 
Milankovitch forcing is nonlinear and that it 
involves some elements of internal forcing. 
However, whether the Milankovitch external 


orbital forcing drives the internal forcing [3], 
phase locks the oscillations of an internally 
driven system [89], or acts as pacemaker for the 
free oscillations of an internally driven system 
[90] is an open question. Furthermore, internal 
forcing mechanisms or feedbacks have 
difficulties in explaining the apparently episodic 
appearance of the 100-kyr cycle [40]. It 
remains to be seen if the wavelet spectra of the 
output data from climate models with internal 
forcing can fit the spectral of the record 
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[72,73]. 

10. PHYSICS BEHIND ICE AGE 
GLACIATION 

The physics behind the marine and 
continental ice age glaciations can be explained 
as follows. The insolation forcing mechanism 
which is capable of reproducing of ice age ice 
sheet is illustrated in Figure 10, in which the 
insolation flux deficit (bfgb'b ) is the physics 
behind the ice age glaciation. The areas ebfe 
and eb'ge represent the regular and prolonged 
glaciation phases respectively; their area ratio is 
about 1:2. It also shows that the minimum 
insolation a b during regular glacial phase 
shifted about 2 Kyr to the minimum insolation 
a'b' during the prolonged glacial phase. While 
the regular glaciation phase lasts only about 4 
Kyr, the prolonged glacial phase lasts longer 
than 8 Kyr. Physically, ice sheet size would 
grow larger if glaciation time is longer and 
glacial phase becomes broader. Therefore, we 
have reason to believe that the negative 
insolation pulse bfgb'b is responsible for the 
100-Kyr cycle of ice age glaciation. Here, we 
identify the insolation flux deficit bfgb'b in 
Figure 10 as the causal mechanism responsible 
for the amplification of an external forcing by 
the orbital resonance dynamics, and arrived in 
this way at a better understanding of the role of 
an externally generated factor in climate 
change. This causal mechanism is appropriate 
and convincing because: 1. The timing of the 
bfgb'b occurrence coincides with the observed 
major climate change events; and 2. The 
insolation flux deficit represented by bfgb'b is 
powerful enough to drive the Earth into an ice 
age. 

The determination of continental 
temperature variations in Siberia by Williams 
et. al. [30] lends fresh support to the theory that 
frequency variation of the obliquity is the main 


mechanism responsible for the late Pleistocene 
climate change. Because Lake Baikal is located 
in the continental interior, its biological 
productivity is sensitive to solar energy 
variations, which are accurately recorded 
through the flux of biogenic silica to the bottom 
sediments. 

The model time series T(t) is shown in 
Fig. 8, superimposed on the time series of the 
oxygen isotope content ($ 8 0) of deep-sea 
sediments and the biogenic silica changes 
(BDP-96-2) of the terrestrial sediments ‘from the 
Baikal Drilling Project (BDP). The T(t) curve 
indeed shows variations of the same order of 
magnitude and phase as seen in the 
paleoclimate records [29, 30, 40]. This implies 
that the marine and the continental 100-Kyr ice 
age cycles are driven by a simple orbital 
mechanism: frequency modulation of the 

obliquity. 

The canonical view of the marine oxygen 
isotope fluctuations is that they are attributed to 
the internal feedback mechanisms in the ocean- 
atmosphere-cryosphere system that can change 
atmosphere temperatures, CO 2 , methane, ice 
sheets, ocean circulation and sea surface 
temperatures [3]. If this view were correct, 
then the far more perplexing mystery may be 
how to explain the continental 100-Kyr ice age 
cycles. Because the continental biological 
productivity is independent of ice sheet size, 
changes in ocean circulation or atmosphere CO 2 
concentration, the physical mechanisms of 
internal forcing for continental climate change 
are far from obvious. The continental 100-Kyr 
puzzle may prove even more mysterious than 
the marine 100-Kyr puzzle [30]. 

11, INTERPRETATION OF GLACIATION 
PUZZLES 

Pulsation of insolation induced by the 
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frequency variations of the Earth's obliquity 
may provide explanations for the 3 glaciation 
puzzles in paleoclimatology. 

The first puzzle in paleoclimatology has 
been understanding why the major glaciation 
cycles of the last 1 million years vary at a 1 00- 
Kyr rhythm. Although the Earth’s orbital 
eccentricity has a 100-Kyr cyclicity, the 
resulting changes in insolation are too small to 
be climatically significant. The second puzzle 
is how to explain the first two million years of 
the northern hemisphere ice ages when global 
ice volume and ocean temperature varied 
predominantly at the 41-Kyr period # of the 
Earth's obliquity. The variation in summer 
insolation at high latitudes is dominated by 
precession. If the Milankovitch view that 
summer radiation at high latitudes exerts the 
ultimate control on ice sheet mass balance is 
correct, then far more precessional variance at 
23 -and 19-Kyr would be observed in high 
latitude climate records of the late Pliocene and 
early Pleistocene. Lastly, we still don't 
understand why the 100-Kyr cycle became 
abruptly dominant about 900 Kyr ago. In the 
past decade, numerous hypotheses have been 
proposed to explain these puzzles; however, 
every explanation has had difficulty accounting 
for the causal mechanisms of orbital or 
stochastical forcing for climate change [3-7, 10, 
11, 13, 14, 16, 17, 91-93]. From a 
climatologist's point of view, the physical 
mechanism of climate change and the rules that 
govern climate transitions between states 
involve ocean-atmosphere-ice processes on 
time scales that are paced by orbital variations. 
Therefore, it is generally believed that only 
through the continued collection and analysis of 
geological data from deep ocean drilling and 
continental drilling can we hope to come to a 
deeper understanding of the mechanisms by 
which the global climate system responds. On 
the contrary, we have shown here that our 


understanding of the dynamics of orbital forcing 
for climate change is far from complete. We 
still have huge gaps in our understanding of the 
first-order orbital resonance effect on climate 
change, gaps which can only be filled by 
fundamental physics. Based on orbital 
resonance theory, this study demonstrates that 
the orbital parameter which paces for the 100- 
Kyr ice age cycle is actually the frequency 
variation of the Earth's obliquity, not orbital 
eccentricity. As shown in Figure 11, 
temperature varies for the past 2 million years 
predomintly at the 4l-Kyr period of the 
obliquity. Also, the modeled 100-Kyr 
temperature variation as shown in Figure 1 1 
becomes dominant abruptly about 900 Kyr ago. 
Why did the 100-Kyr cycle only appear about 
0.9 million years ago? Because since that time 
cooperative phenomena between the frequency 
and amplitude variations of the obliquity have 
induced a resonance effect on climate response 
to insolation, so as to shift the 41-Kyr cycle into 
the 100-KyT cycle. 

12. CONCLUDING REMARKS 

The study of how frequency variation of 
the Earth's obliquity came to be related with the 
paleoclimatic changes has recently moved into 
a new phase of scientific inquiry and discovery. 
Notable advances in atmospheric science have 
given rise to new ways of insolation calculation 
to solve the perplexing and enduring glaciation 
puzzles in paleoclimatology. The focus of this 
review is on the problem in insolation theory of 
climate change. The problem is that all ice 
models were introduced in the climate system 
without ' the appropriate time-dependent 
response to the energetics of the climate 
system. In order to make up for this deficiency 
in the climate system, we have introduced a 
pulsation term of insolation as a device for 
illustrating the physics behind the glacial cycles. 
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This has been done in terms of a bipolar pulse 
modulation train c(t) of insolation. According 
to theories of resonance, insolation pulsation 
has been applied to characterize the time- 
dependent response of ice models in the climate 
system. Incorporating the amplitude-frequency 
resonance effect of the obliquity, climate model 
simulations have delivered an ice age 
chronology that is in close accord with <5> 8 0 and 
BDP proxy data. This methodology of 
insolation pulsation is developed from a 
convolution theory of the oon-stationary time 
series of insolation by the Hibert transformation 
in order to find out the three glaciation puzzles 
in paleoclimatology. The success of the model 
simulation results is remarkable because they 
can provide explanations to account for the 
origin, structure and spectrum of the observed 
major cl|mate changes during the past 2 million 
years. Ik conclusion, we claim that frequency 
modulation effect of the Earth's obliquity on the 
incoming solar radiation may be the key to 
solving the century-long glaciation puzzles in 
climatology. 
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POPULAR SUMMARY 


Astonomical Background: 


* Frequency modulation of the Earth's obliqity is the 
major cause of climate change. It was misleading to call 
upon other astronomical forcing such as orbital eccentricity, 
inclination and precession to explain the origin of the 
ice-age cycle which is a mystery in science. 

Geophysical Speculations : 

/ 

* The varying frequency of the obliquity can induce pulses 
in the incoming solar radiation. Glacial and interglacial 
climate events are the consequence of insolation pulsation. 

It is unnecessay to speculate internal feedbakes such as 
ocean circulations and atmospheric CO 2 concentrations in 
the Earth system to understand the glaciation puzzles. 

Geological Data: 

* Results from model simulations are in good agreement with 
geological climate records for the past 2 million years . 


Climate Prediction: 


* Extension of model simulation reveals that the current 
warning trend of the climate is almost over and will give 
way to forecast a small ice age. 



